In [3]:

    

# plotting
import holoviews as hv
%load_ext holoviews.ipython
hv.notebook_extension('bokeh')

# scientific python usage
import numpy as np
from   numpy import exp,cos,sin,pi,tan,sqrt
from scipy.stats import norm
from __future__ import division

# dedist package
from DeDist import dedist


    

    




The holoviews.ipython extension is already loaded. To reload it, use:
  %reload_ext holoviews.ipython






    












    















    















  






  
  

In [4]:

    

def fun(x,x_,par):
    ''' Gaussian tuning curve
    
    Parameters
    ----------
    x,x_ : arrays
        Preferred and actual stimuli
    par : array or list
        Tuning curve parameters, should be as [a,s], with 
        a being the max response, and s**2 the variance.
    
    Returns
    -------
    array
         of form a*exp(-(x-x_)**2/(2*s**2))
    '''
    return par[0]*exp(-(x-x_)**2/(2*par[1]**2))


    

In [5]:

    

# number of neurons
N = 64

# set of orientations
A = np.linspace(-pi,pi,N+1)[:-1]

# gaussian noise covariance matrix
c = 1
ro = 1
sigma = 0.5
cov = sigma**2*c*exp(-abs(A-A[:,None])/ro)
cov[range(len(A)), range(len(A))] = np.ones(len(A))*sigma**2

# Stimulus sampling (at what stimuli the probability dist should be calculated)
# This needs to be at least as high as N. The stimulus sampling in the case of 
# mirrored functions (such a function dependent on the absolute value of the stimulus)
# the stimulus samples should only include half of the possible inputs
ns = 64
As = np.linspace(-pi,pi,ns+1)[:-1]

# set parameters
par = [1,1]


    

In [15]:

    

# number of realizations
reals = 1000

# set stimulus index (such that the stimulus value is As[stim])
stim = 32

# define noise
n = np.random.multivariate_normal(np.zeros(N), cov, size=(reals))

# find all possible responses across sample stimuli
F = fun(A,As[:,None],par)

# find real response and add noise
f = F[stim,:]
r = f+n

# find total squared errors across possible stimuli
F_errors = r - F[:, None]
tmp = np.tensordot(F_errors, np.linalg.inv(cov), axes=1)
Errors = np.einsum('ijk,ijk->ij', tmp, F_errors)

# find the minimum error stimuli for each noise realization
sol = As[Errors.argmin(axis=0)]

# plot decoding distribution histogram
dA = (As[1]-As[0])/2
freq, edges = np.histogram(sol, bins=As-dA)
hv.Histogram(freq/freq.sum(),edges)


    

    
Out[15]:






    

        

    

In [16]:

    

p, means, covs = dedist.est_p_cor(fun,As[stim],par,cov,A,As,verbose=False, full_return=True)
hv.Histogram(freq/freq.sum(),edges)*hv.Curve(zip(As,p))


    

    
Out[16]: